clear
global nbar Nmean N2mean Nmax X Y R2 prXY XY XY4 dr

% optimalizace kaskad

% Definice funkci
%%%%%%%%%%%%%%%%%%%%%%
% termalni vstup, hleda se chi a alpha aby to melo nejvetsi nmean-Delta n
function y=Dolik(x)
 global nbar
 alpha=x(1);
 c=cos(alpha);
 s=sin(alpha);
 chi=x(2);
 G=c^2*s^2*sin(chi)/nbar^2/((1/nbar+2*s^2*sin(chi/2)^2)^2+s^4*sin(chi)^2)^2;
 nbarnew=c^2*nbar+G;
 zavorka1=(1/nbar+2*s^2*sin(chi/2)^2)^2+s^4*sin(chi)^2;
 G1=c^2*s^2*(1-2*c^2*nbar)*sin(chi)/nbar^2/zavorka1^2;
 G2=-c^4*s^4*sin(chi)^2/nbar^4/zavorka1^4;
 G3=4*c^4*s^2*sin(chi)*(1/nbar+2*s^2*sin(chi/2)^2)/nbar^2/zavorka1^3;
 zavorka2=(1/nbar+2*s^2*sin(chi)^2)^2+s^4*sin(2*chi)^2;
 G4=c^4*s^4*sin(2*chi)^2/nbar^2/zavorka2^3;
 Deltan=sqrt(c^4*nbar^2+c^2*nbar+G1+G2+G3+G4);
 y = -(nbarnew-Deltan);
% y=-nbarnew;
endfunction

% pro nalezeni bet a gam v MaxEnt funkci
function y=OdchylEnt(x)
 global Nmean N2mean Nmax
 bet=x(1);
 gam=x(2);
 prob=exp(-([0:Nmax]-bet).^2/2/gam^2);
 prob=prob/sum(prob);
 nmeanp=prob*[0:Nmax]';
 n2meanp=prob*([0:Nmax].^2)';
 y=(nmeanp-Nmean)^2+(n2meanp-N2mean)^2;
endfunction






nbar=100;
%Dolik([0.2;0.3])

pomodhad=[0.3;0.1];
vals=fminunc('Dolik',pomodhad)

% optimalizovane hodnoty
alpha=vals(1);
chi=vals(2);
c=cos(alpha);
s=sin(alpha);
G=c^2*s^2*sin(chi)/nbar^2/((1/nbar+2*s^2*sin(chi/2)^2)^2+s^4*sin(chi)^2)^2;
nbarnew=c^2*nbar+G;
zavorka1=(1/nbar+2*s^2*sin(chi/2)^2)^2+s^4*sin(chi)^2;
G1=c^2*s^2*(1-2*c^2*nbar)*sin(chi)/nbar^2/zavorka1^2;
G2=-c^4*s^4*sin(chi)^2/nbar^4/zavorka1^4;
G3=4*c^4*s^2*sin(chi)*(1/nbar+2*s^2*sin(chi/2)^2)/nbar^2/zavorka1^3;
zavorka2=(1/nbar+2*s^2*sin(chi)^2)^2+s^4*sin(2*chi)^2;
G4=c^4*s^4*sin(2*chi)^2/nbar^2/zavorka2^3;
Deltan=sqrt(c^4*nbar^2+c^2*nbar+G1+G2+G3+G4);

nbarnew-Deltan;

%Aproximace fotodistribucni funkce pomoci MaxEnt
Nmax=1200;
Ns=0:Nmax;

Nmean=nbarnew;
N2mean=Nmean^2+Deltan^2;

%Nmean=118.14;
%Deltan=91.497;
%Nmean=127.12;
%Deltan=95.556;
%Nmean=140.22  
%Deltan=103.72;
%Nmean=156.76;
%Deltan=113.90;
%Nmean=177.74;
%Deltan=127.28;
%Nmean=204.17;
%Deltan=143.38;
Nmean= 238.78
Deltan= 165.43
N2mean=Nmean^2+Deltan^2;

%pomodhad=[-10;100];
pomodhad=[100;250];
vals=fminunc('OdchylEnt',pomodhad);
OdchylEnt(vals)
pomodhad=vals;
vals=fminunc('OdchylEnt',pomodhad);
OdchylEnt(vals)
pomodhad=vals;
vals=fminunc('OdchylEnt',pomodhad);
OdchylEnt(vals)

bet=vals(1)
gam=vals(2)
probMEnt=exp(-([0:Nmax]-bet).^2/2/gam^2);
probMEnt=probMEnt/sum(probMEnt);
nmeanp=probMEnt*[0:Nmax]'
n2meanp=probMEnt*([0:Nmax].^2)';
DeltanNew=sqrt(n2meanp-nmeanp^2)


figure(1)
plot(Ns,probMEnt,'r')



% Calculation of new mean photon number and fluctuation out of the Glauber Sudarshan
%q1=-3.0880e-03;
%q2=  -1.0211e-03;
%q3=  -2.3559e-04;
%q4=  -8.0570e-06;
%q1=-2.3663e-03;
%q2=  -1.0011e-03;
%q3=   2.0870e-05;
%q4=  -2.0653e-05;
%q1=-2.8251e-03;
%q2=   9.0640e-04;
%q3=  -8.4936e-05;
%q4=  -1.7042e-05;
%q1=-3.5010e-03;
%q2=   2.1207e-03;
%q3=  -1.0831e-04;
%q4=  -1.5621e-05;
%q1=-3.9547e-03;
%q2=   3.0690e-03;
%q3=  -1.5025e-04;
%q4=  -1.2327e-05;
%q1=-5.4089e-03;
%q2=   2.8786e-03;
%q3=  -5.0374e-05;
%q4=  -1.3404e-05;
q1=-6.0062e-03;
q2=   3.3154e-03;
q3=  -7.6323e-05;
q4=  -1.0003e-05;


rmax=sqrt(Nmax);
dr=rmax/200;
r=0:dr:rmax;
pr=exp(q4*r.^4+q3*r.^3+q2*r.^2+q1*r).*r;
pr=pr/sum(pr)/dr;

alpha=0.15
chi=0.15
s=sin(alpha)
c=cos(alpha)

[X,Y]=meshgrid(r,r);
R2=X.^2+Y.^2;
XY=X.*Y;
fcePom=exp(-(2*s^2*R2*sin(chi/2)^2)).*(R2+2*XY.*besselj(1,2*s^2*XY*sin(chi)));
prXY=pr'*pr;

nbarNew=c^2/2*sum(sum(prXY.*fcePom))*dr^2
normal=sum(sum(prXY))*dr^2

% vypocet variance n
XY4=X.^4+Y.^4+4*X.^2.*Y.^2;
pomTerm1=4*XY.*R2.*exp(-(2*s^2*R2*sin(chi/2)^2)).*besselj(1,2*s^2*XY*sin(chi));
pomTerm2=XY.^2.*exp(-(2*s^2*R2*sin(chi)^2)).*besselj(2,2*s^2*XY*sin(2*chi));
momAs=c^4/4*sum(sum(prXY.*(XY4+pomTerm1+pomTerm2)))*dr^2;

DeltanNew=sqrt(momAs+nbarNew-nbarNew^2)

figure(2)
mesh(r,r,prXY.*(XY4+pomTerm1+pomTerm2))


% optimalizace nbar - Deltan
function y=Dolik2(x)
 global X Y R2 prXY XY XY4 dr
 alpha=x(1);
 c=cos(alpha);
 s=sin(alpha);
 chi=x(2);
 fcePom=exp(-(2*s^2*R2*sin(chi/2)^2)).*(R2+2*XY.*besselj(1,2*s^2*XY*sin(chi)));
 nbarnew=c^2/2*sum(sum(prXY.*fcePom))*dr^2;
 pomTerm1=4*XY.*R2.*exp(-(2*s^2*R2*sin(chi/2)^2)).*besselj(1,2*s^2*XY*sin(chi));
 pomTerm2=XY.^2.*exp(-(2*s^2*R2*sin(chi)^2)).*besselj(2,2*s^2*XY*sin(2*chi));
 momAs=c^4/4*sum(sum(prXY.*(XY4+pomTerm1+pomTerm2)))*dr^2;
 Deltan=sqrt(momAs+nbarnew-nbarnew^2);
% y = -(nbarnew-0.9*Deltan);
 y = -(nbarnew-Deltan);
% y=-nbarnew;
endfunction

pomodhad=[0.47;0.03];
pomodhad=[0.2;0.2];
pomodhad=[0.4;0.02];
vals=fminunc('Dolik2',pomodhad)


 alpha=vals(1);
 c=cos(alpha);
 s=sin(alpha);
 chi=vals(2);
 fcePom=exp(-(2*s^2*R2*sin(chi/2)^2)).*(R2+2*XY.*besselj(1,2*s^2*XY*sin(chi)));
 nbarnew=c^2/2*sum(sum(prXY.*fcePom))*dr^2
 pomTerm1=4*XY.*R2.*exp(-(2*s^2*R2*sin(chi/2)^2)).*besselj(1,2*s^2*XY*sin(chi));
 pomTerm2=XY.^2.*exp(-(2*s^2*R2*sin(chi)^2)).*besselj(2,2*s^2*XY*sin(2*chi));
 momAs=c^4/4*sum(sum(prXY.*(XY4+pomTerm1+pomTerm2)))*dr^2;
 Deltan=sqrt(momAs+nbarnew-nbarnew^2)

